* This do-file, run by ssdi_master.do, generates several data quality checks
*	as well as the final summary stats table. 

* NOTE: this do-file assumes use of SSB v.6.02

******* Section 1: Preliminaries *******

set more off // other base settings in master do-file

*suppress output for most data cleaning

qui{

******* Section 2: Creating the Sample *******

*** Data Restriction 1: restrict to the relevant SIPP panel years

qui keep if inlist(panel, 1996, 2001, 2004, 2008)
gen year = panel

/**** Make family linkage to identify which households have dependents
 (and are thus receiving dependent benefits). */

qui preserve

qui keep person_id spouse_person_id

qui ren (person_id spouse_person_id) (id_p id_s)

qui gen fid = _n

qui reshape long id_, i(fid) j(person) str

qui collapse (firstnm) fid, by(id)

qui drop if mi(id_)

qui ren id_ person_id
qui ren fid fam_id

qui save famids.dta, replace

qui restore

*merge the family ids onto the full dataset

qui merge 1:1 person_id using famids.dta
qui drop _merge

**calculate if there are any dependents in the family for female spouses

*NOTE: last_admin_birthdate = birthdate of youngest child if observation female, 
*	15-65 during fertility history module (which was asked in all panels)

gen youngest_age = (sipp_panel_beg_date-last_admin_birthdate)/365.25 // will be missing if male or didn't have a child, age in years otherwise
gen hasdepkids = 1 if youngest_age<18 // will be missing if male, didn't have a child, or child was over 18; 1 otherwise

*apply this determination to both males and females in the family
bys fam_id: egen hasdep = min(hasdepkids) // if female spouse has hasdepkids==1, then will set this for male spouse as well

/*** Data quality change: 
per email from Jordan at Census, need to adjust the ssdi amounts from the 
MBR file. This will be fixed after SSB v.6.02. */	

replace mbr_ssdi_benefit_totamt_1 = mbr_ssdi_benefit_totamt_1/10

*** generate age and panel year variables 

*age
gen sipp_beg_age = (sipp_panel_beg_date - birthdate)/365.25 
gen sipp_end_age = (sipp_panel_end_date - birthdate)/365.25

*panel years
gen panel_startyr = year(sipp_panel_beg_date)
gen panel_finyr = year(sipp_panel_end_date)

/*** Restrictions from the paper, following the admin data sample construction ****
NOTE: In the core SSA data we enforce four restrictions, but we can only enforce two of them on 
the SSB data. */

/*Data Restriction 2: disability case not reduced for age
NOTE: We want to drop this group altogether from the sample
NOTE2: each of benefit values are indicators for receiving benefits conditional on applying, missing o/w */

gen rest1= (inlist(1, mbr_retire_benefit, mbr_agedsp_benefit, mbr_widowsp_benefit) & sipp_end_age < 60) // indicator if receiving retirement/survivors insurance while under age 60
la var rest1 "if receiving retirement/survivors insurance while under age 60"
drop if rest1==1

/*Data Restriction 3: there is a single date of filing
NOTE1: Since we do not observe date of filing, we instead restrict to cases 
	where there's only one occurrence of eligibility. 
NOTE2: If never applied, then kept in the "non-SSDI" group; if applied multiple times
	or applied out of order, then drop from the sample altogether. For the last two cases, 
	I create a flag here; I summarize and drop immediately before the summary stats 
	so we will know the number of observations as a data quality measure. */

* generate multiple app flag (=1 if observations have mbr_ssdi_entitled_`i' = 1 for i=2, 3, or 4 AND mbr_ssdi_entitled_1 = 1)

* SSB variable is missing if didn't apply for benefits, 0 if denied, 1 if accepted; only care about last case, so generate an indicator for whether accepted
forv i=1/4 {
recode mbr_ssdi_entitled_`i' (0 . = 0), gen(mbr_ssdi_entitled_ind_`i') // indicator for whether application i yielded entitlement to benefits
}

* now generate flag for multiple eligible applications

egen mult_elig_apps = rowtotal(mbr_ssdi_entitled_ind_1 mbr_ssdi_entitled_ind_2 mbr_ssdi_entitled_ind_3 mbr_ssdi_entitled_ind_4) // counts the number of eligible apps for i=1/4; never missing
gen mult_flag = (mult_elig_apps>1) // only >1 if (at least) two of mbr_ssdi_entitled_ind_i is ==1
la var mult_flag "=1 if eligible application after the first"

* generate flag for out of sequence observations

gen oos_flag = (mbr_ssdi_entitled_ind_1 ==0 & inlist(1, mbr_ssdi_entitled_ind_2, mbr_ssdi_entitled_ind_3, mbr_ssdi_entitled_ind_4)) // flag =1 if first app not elig but one of the others is
la var oos_flag "=1 if first not eligible but other applications are"

*later we'll summarize the dummies, and want to know their proportions as a fraction of the SSDI population.

/* Identify SSI recipients
NOTE: ssr_ssi_benefit_fed_totamt = missing if no documented SSI receipt from federal sources 
	(check federal sources only in order to be consistent with the admin data) */
	
/* identify individuals that received SSI before/during the panel (these are "SSI recipients"
for our purposes, as they will have SSI affect their SIPP characteristics); o/w not SSI. 
NOTE1: people who received SSI before the panel should be classified as SSI recipients (serial beneficiaries) */

gen ssi=.
replace ssi = 0 if mi(ssr_ssi_benefit_fed_totamt) // ie if never received benefits
replace ssi = 1 if (ssr_ssi_benefit_fed_totamt>0 & ~mi(ssr_ssi_benefit_fed_totamt)) // had received SSI at any point
replace ssi = 0 if (ssr_ssi_first_pmt_dt > sipp_panel_end_date) // was recipient after the panel ended, so okay to include in the sample. NOTE: if mi(ssr_ssi_first_pmt_dt), then will set ssi==0; but this group already has no SSI benefits so should be no change

/* identify PIA variables
As soon as we restrict to people with one date of eligibility, all people have 
only one PIA value, which is stored in mbr_ssdi_benefit_totamt_1 (eligible), 
phus_ssdi_benefit_totamt_1 (paid). Moreover, the eligibility determination date 
is currently stored in mbr_ssdi_doed_1. */

replace mbr_ssdi_doed_1 = . if mbr_ssdi_entitled_1==0 // set to no eligibility date if denied for benefits
replace mbr_ssdi_doed_1 = . if mbr_ssdi_benefit_totamt_1==0 // set to no eligibility date if you have no recorded benefits (small amount of cases)

*generate dates of eligibility 

form mbr_ssdi_doed_1 %td
gen pia_year = year(mbr_ssdi_doed_1)
gen pia_month = month(mbr_ssdi_doed_1)
gen pia_yyyymm = pia_year*100 + pia_month

gen pia_nom = mbr_ssdi_benefit_totamt_1 // Take PIA to be the eligible amount, per our codebook

*store MBR SSDI diagnosis groups

gen piareason = mbr_ssdi_dig_group_1
la var piareason "MBR SSDI: Diagnosis group"
la def lmbr_ssdi_dig_group_1 ///
	1  "Diseases of the blood" ///
	2  "Circulatory system" ///
	3  "Congenital anomalies" ///
	4  "Digestive system" ///
	5  "Endocrine, nutritional, and metabolic diseases" ///
	6  "Genitourinary system" ///
	7  "Infectious and parasitic diseases" ///
	8  "Injuries" ///
	9  "Mental disorders" ///
	10  "Mental retardation" ///
	11  "Musculoskeletal system" ///
	12  "Neoplasms" ///
	13  "Nervous system and sense organs" ///
	14  "Respiratory system" ///
	15  "Skin and subcutaneous tissue" ///
	16  "Other or unknown"
la val piareason lmbr_ssdi_dig_group_1

*NOTE: can't replicate two parts of the code for admin data-- variables for date 
*  of application for SSDI, date of filing for SSDI, both of which we don't observe. 

* generate variable for "age when application was filed" (use age when declared 
*	eligible for benefits instead, per our codebook)

gen age_elig = int((mbr_ssdi_doed_1-birthdate)/365.25) // missing only if never applied
gen age_filing = age_elig
la var age_elig "Age declared eligible for SSDI"

* race identifiers (1 = white, 2 = black, 3 = other)

recode race (2 3 = 0), gen(white)
recode race (1 3 = 0) (2 = 1), gen(black)
recode race (1 2 = 0) (3 = 1), gen(otherrace)

*identify the IME year from the PIA yyyymm

gen imeyear=.

forv i=1992/2012 { //this is the lower bound year and the imeyear
loc j = `i' + 1 // this is the upper bound year

replace imeyear=`i' if `i'06<=pia_yyyymm & pia_yyyymm<=`j'05 // observation is in a year `i' if between June of that year, May of next
}

*create index for the national wage index and CPI. I'll "merge" using locals and 
*	a loop instead of adding an extra dataset.

*make empty variables to replace into
gen awi_year=.
gen wage_index=.
gen cola=.
gen cola_cum=.

*now make arrays, index for them
loc imeyear_array 1992 1993 1994 1995 1996 1997 1998 1999 2000 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012
loc awi_year_array 1990 1991 1992 1993 1994 1995 1996 1997 1998 1999 2000 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010
loc wage_index_array 21027.98 21811.6 22935.42 23132.67 23753.53 24705.66 25913.9 27426 28861.44 30469.84 32154.82 32921.92 33252.09 34064.95 35648.55 36952.94 38651.41 40405.48 41334.97 40711.61 41673.83
loc cola_array 3 2.6 2.8 2.6 2.9 2.1 1.3 2.5 3.5 2.6 1.4 2.1 2.7 4.1 3.3 2.3 5.8 0 0 3.6 1.7
loc cola_cum_array 1.683504961 1.634470836 1.593051497 1.549660989 1.510390828 1.467823934 1.437633628 1.419184233 1.384569983 1.337748776 1.303848709 1.285846853 1.259399465 1.226289644 1.177991973 1.14036009 1.114721496 1.053612 1.053612 1.053612 1.017
loc yearslength: word count `imeyear_array'

forv i=1/`yearslength' { //think of this as picking out the 1x5 data vector for each year
*store each number in seperate locals
loc imeyear: word `i' of `imeyear_array'
loc awi_year: word `i' of `awi_year_array'
loc wage_index: word `i' of `wage_index_array'
loc cola: word `i' of `cola_array'
loc cola_cum: word `i' of `cola_cum_array'

replace awi_year=`awi_year' if imeyear==`imeyear'
replace wage_index=`wage_index' if imeyear==`imeyear'
replace cola=`cola' if imeyear==`imeyear'
replace cola_cum=`cola_cum' if imeyear==`imeyear'
}

/*Deflating the PIA on the later six months based only on the CPI increase for that year, 
because PIA is increased even though AIME stays the same*/

gen pia_step = .
replace pia_step = pia_nom if 6<=pia_month & pia_month<=11
replace pia_step = pia_nom/(1+(cola/100)) if ~(6<=pia_month & pia_month<=11)

* Converting to 2013 CPI values
gen pia = pia_step*cola_cum 

* Add variables for when first payments are made
gen startpay_year = year(phus_ssdi_benefit_stdate_1)
gen startpay_month = month(phus_ssdi_benefit_stdate_1)
gen startpay_yyyymm = startpay_year*100 + startpay_month
gen startpay = phus_ssdi_benefit_stdate_1

/* Data Restriction 4: age restriction should be: 21-61 throughout the SIPP panel AND 
	age at eligibility over 21 (for those who were eligible)
NOTE: age_elig is missing if individual never eligible for SSDI. */

keep if sipp_beg_age >=21 & sipp_end_age<=61 // keep if people are aged 21-61
drop if age_elig<21 & ~mi(age_elig) // for those who were eligible, drop if less than 21 at elig

/* Data Restriction 5: length of time before payments
*Restrict analysis sample to those with date of eligibility
	being less than 65 months from the date payments start. 
*NOTE: 365.25 days/12 months = 30.4375 days in a month on average */

gen months_payment_applic = (phus_ssdi_benefit_stdate_1-mbr_ssdi_doed_1)/30.4375 // months in between
drop if months_payment_applic > 4 & ~mi(months_payment_applic)

* Data Restriction 6: need to be on DI before panel starts
drop if (phus_ssdi_benefit_stdate_1 > sipp_panel_beg_date) & ~mi(phus_ssdi_benefit_stdate_1)

* Calculate how many months between death and first payment
gen death_months = (deathdate-startpay)/(365.25/12)

/***** identify bend points using the PIA distribution from the admin data
For each panel, we have bend points at the following points in the cdf:
	LBP: 0.0055-0.1184
	FMBP: 0.1294-0.4756
	UBP: 0.7658-0.9156 
*/

sort pia

qui levelsof year // handy way to grab the four years
foreach yr in `r(levels)' { 
cumul pia if year==`yr' & ~mi(pia) & ssi==0, g(dist_`yr') eq // estimate empirical cdf of PIA for each panel, store in dist_`yr'

*now identify the bend points
gen l_`yr' = 1 if  0.0055<= dist_`yr' & dist_`yr'<=0.1184
gen u_`yr' = 1 if  0.7658<= dist_`yr' & dist_`yr'<=0.9156
gen f_`yr' = 1 if  0.1294<= dist_`yr' & dist_`yr'<=0.4756 & hasdep==1
}

*now sum across panels to generate the sample

gen group_bp1 = 1 if inlist(1, l_1996, l_2001, l_2004, l_2008) // =1 if in the LBP part of the distribution
gen group_bp2 = 1 if inlist(1, u_1996, u_2001, u_2004, u_2008)
gen group_bpf = 1 if inlist(1, f_1996, f_2001, f_2004, f_2008)

/* Make 1-missing groups for the summary stats table. We want the following groups in the table:
1. DI bend points (already made)
2. all SSDI recipients ("all_ssdi")
3. all non-SSDI recipients ("all_no_ssdi")
4. SSDI recipients, no SSI ("ssdi_no_ssi")
5. SSI recipients, no DI ("ssi_no_ssdi")
6. DI and SSI joint recipients ("ssi_and_ssdi")
7. no SSI, no DI ("no_ssi_no_ssdi")
8. all SSI recipients ("all_ssi")
*/

gen all_ssdi = 1 if (pia>0 & ~mi(pia)) // 	 all DI recipients
gen all_no_ssdi = 1 if mi(all_ssdi) // 				 all non-DI recipients
gen ssdi_no_ssi = 1 if (all_ssdi==1 & ssi==0) // 	 SSDI recipients, no SSI
gen ssi_no_ssdi = 1 if (mi(all_ssdi) & ssi==1) //	 SSI recipients, no SSDI
gen ssi_and_ssdi = 1 if (all_ssdi==1 & ssi==1) //	 SSI and SSDI recipients
gen no_ssi_no_ssdi = 1 if (mi(all_ssdi) & ssi==0) // no SSI, no SSDI
gen all_ssi = 1 if ssi==1

******* Section 3: Prep for Summary Stats*******

*demographic variables

qui tab state, gen(state_)
ren hispanic hisp
recode educ_5cat (2 3 4 5 = 0), gen(nohs)
recode educ_5cat (1 3 4 5 = 0) (2 = 1), gen(hsdeg)
recode educ_5cat (1 2 4 5 = 0) (3 = 1), gen(somecoll)
recode educ_5cat (1 2 3 5 = 0) (4 = 1), gen(colldeg)
recode educ_5cat (1 2 3 4 = 0) (5 = 1), gen(graddeg)

*enter CPI-U-RS into the table

gen cpi_urs_2013 = .

*I have the full series in case we need to adjust later 
loc cpi_yrlist 1996 2001 2004 2008 
loc cpi_vals 0.675920514 0.759789597 0.890122735 

loc cpi_totyrs: word count `cpi_yrlist'

forv i=1/`cpi_totyrs'{
loc cpiyr: word `i' of `cpi_yrlist'
loc cpival: word `i' of `cpi_vals'

replace cpi_urs_2013 = `cpival' if year == `cpiyr'
}

*generate indicators for the different years
foreach yr in 1996 2001 2004 2008 {
	gen yr`yr' = 1 if year == `yr'
}

*for all variables observed in a series, generate variable that holds value for 
*  the first variable observed in the panel year (eg for yearly, panel year 
*  variable; for monthly, first variable in panel year).

gen totinc = totinc_1
gen totearn = totearn_1
gen fsamt = fsamt_1
gen hicov = hicov_1
gen hiemp = hiemp_1
gen afdcamt = afdcamt_1
gen vcamt = vcamt_1
gen wcamt = wcamt_1

gen totearn_ser =.
gen total_der_fica =.
gen total_der_nonfica =.

qui levelsof panel

*the addition of 1.1081 is to scale it from 2013 to 2020 dollar values
foreach var in totearn_ser total_der_fica total_der_nonfica {
	foreach year in `r(levels)' {
		replace `var' = `var'_`year' * 1.1081 if panel==`year'
	}
}

local lhs own_home db_pension dc_pension pension_in_scope_empl life_ins_1 hicov hiemp

loc nom_lhs nonhouswealth homeequity totnetworth totinc totearn totearn_ser total_der_fica total_der_nonfica fsamt afdcamt vcamt wcamt

*Adjust to 2020 dollars using cpi_urs
foreach var in `nom_lhs' {
	capture drop real_`var'
	gen r_`var' = (`var'/cpi_urs) 
	_crcslbl r_`var' `var'
}

loc real_lhs r_nonhouswealth r_homeequity r_totnetworth r_totinc r_totearn r_totearn_ser r_total_der_fica r_total_der_nonfica r_fsamt r_afdcamt

}

log using "gs_${gsnum}${impflag}.txt", t replace

** data quality checks

*mult_flag, oos_flag

*percentage of total pop
qui su mult_flag
di "mean: `r(mean)'; N: `r(N)'"
qui su oos_flag
di "mean: `r(mean)'; N: `r(N)'"

*percentage of ssdi pop
qui su mult_flag if ssdi_no_ssi==1
di "mean: `r(mean)'; N: `r(N)'"
qui su oos_flag if ssdi_no_ssi==1
di "mean: `r(mean)'; N: `r(N)'"

*percentage at each of the BP's
qui su mult_flag if group_bp1==1
di "mean: `r(mean)'; N: `r(N)'"
qui su oos_flag if group_bp1==1
di "mean: `r(mean)'; N: `r(N)'"

qui su mult_flag if group_bp2==1
di "mean: `r(mean)'; N: `r(N)'"
qui su oos_flag if group_bp2==1
di "mean: `r(mean)'; N: `r(N)'"

qui su mult_flag if group_bpf==1
di "mean: `r(mean)'; N: `r(N)'"
qui su oos_flag if group_bpf==1
di "mean: `r(mean)'; N: `r(N)'"

*drop the flagged obs
drop if inlist(1, oos_flag)
drop if inlist(1, mult_flag)

** pia means for each of the groups

qui su mbr_ssdi_benefit_totamt_1 if group_bp1==1
di "mean: `r(mean)'; N: `r(N)'"
qui su mbr_ssdi_benefit_totamt_1 if group_bp2==1
di "mean: `r(mean)'; N: `r(N)'"
qui su mbr_ssdi_benefit_totamt_1 if group_bpf==1
di "mean: `r(mean)'; N: `r(N)'"
qui su mbr_ssdi_benefit_totamt_1 if all_ssdi==1
di "mean: `r(mean)'; N: `r(N)'"
qui su mbr_ssdi_benefit_totamt_1 if all_no_ssdi==1
di "mean: `r(mean)'; N: `r(N)'"
qui su mbr_ssdi_benefit_totamt_1 if ssdi_no_ssi==1
di "mean: `r(mean)'; N: `r(N)'"
qui su mbr_ssdi_benefit_totamt_1 if ssi_no_ssdi==1
di "mean: `r(mean)'; N: `r(N)'"
qui su mbr_ssdi_benefit_totamt_1 if ssi_and_ssdi==1
di "mean: `r(mean)'; N: `r(N)'"
qui su mbr_ssdi_benefit_totamt_1 if no_ssi_no_ssdi==1
di "mean: `r(mean)'; N: `r(N)'"
qui su mbr_ssdi_benefit_totamt_1 if all_ssi==1
di "mean: `r(mean)'; N: `r(N)'"

qui su pia if group_bp1==1
di "mean: `r(mean)'; N: `r(N)'"
qui su pia if group_bp2==1
di "mean: `r(mean)'; N: `r(N)'"
qui su pia if group_bpf==1
di "mean: `r(mean)'; N: `r(N)'"
qui su pia if all_ssdi==1
di "mean: `r(mean)'; N: `r(N)'"
qui su pia if all_no_ssdi==1
di "mean: `r(mean)'; N: `r(N)'"
qui su pia if ssdi_no_ssi==1
di "mean: `r(mean)'; N: `r(N)'"
qui su pia if ssi_no_ssdi==1
di "mean: `r(mean)'; N: `r(N)'"
qui su pia if ssi_and_ssdi==1
di "mean: `r(mean)'; N: `r(N)'"
qui su pia if no_ssi_no_ssdi==1
di "mean: `r(mean)'; N: `r(N)'"
qui su pia if all_ssi==1
di "mean: `r(mean)'; N: `r(N)'"

log close

*back to suppressing output
qui{

/* Now it's time to calculate summary stats by group and by year 
For example, the average total expenditures in the last quarter in 1999 just among households within the bandwidth at the lower bend point in 1999. 
We save each of those statistics in a local and pull them later to put into an Excel sheet. */

*generate some more variables for the table:

gen married = ~mi(spouse_person_id)
la var married "=1 if married"
gen nkids = own_kids_ever
la var nkids "number of kids"
recode piareason (9 10 = 1) (1 2 3 4 5 6 7 8 11 12 13 14 15 16 = 0), gen(reason_mental)
la var reason_mental "=1 if mental SSDI diagnosis group"
recode piareason (11 = 1) (1 2 3 4 5 6 7 8 9 10 12 13 14 15 16 = 0), gen(reason_musculo)
la var reason_musculo "=1 if musculo SSDI diagnosis group"
recode piareason (1 2 3 4 5 6 7 8 12 13 14 15 16 = 1) (9 10 11 = 0), gen(reason_other)
la var reason_other "=1 if other SSDI diagnosis group"
gen fs_ind = (fsamt>0) & ~mi(fsamt)
la var fs_ind "=1 if received food stamps"

*fraction with no earnings
gen no_earn = (inlist(totearn, 0, .)) // note that totearn is never missing (imputed if missing)
la var no_earn "=1 if no earnings"

/* create variable for fraction of DI subgroup with no non-DI income, calculated by checking if any source of income in the dataset is positive
Sources of income available in the SSB: earnings (totearn), food stamps (fsamt), AFDC assistance (afdcamt), veterans compensation (vcamt), 
	SSI (will use prior SSI determination), workers compensation (wcamt) */

loc inc_types totearn fsamt afdcamt vcamt wcamt

* convert amounts to indicators so we can use inlist
foreach var in `inc_types' {
gen `var'_ind = .
replace `var'_ind = 1 if `var'>0 & ~mi(`var')
}

gen inc_allDI = all_ssdi // restricts to just the DI population
replace inc_allDI = 0 if inlist(1, totearn_ind, fsamt_ind, afdcamt_ind, vcamt_ind, wcamt_ind, ssi) & ~mi(all_ssdi) // for DI recips, set to zero if also getting other income
la var inc_allDI "=1 if on DI and no non-DI income"

*variable labels for the summary stats
la var male "Male"
la var sipp_beg_age "Age at panel start"
la var white "White"
la var black "Black"
la var hisp "Hispanic"
la var otherrace "Other Race"
la var nohs "No HS deg."
la var hsdeg "HS deg."
la var somecoll "Attended college, no deg."
la var colldeg "College deg."
la var graddeg "Graduate deg."
la var r_nonhouswealth "Non-housing wealth"
la var own_home "Owns home"
la var r_homeequity "Home Equity"
la var r_totnetworth "Total Net Worth"
la var r_totinc "Total Income, First Month"
la var r_totearn "Total Earnings, First Month"
la var r_totearn_ser "Total Yearly Earnings, First Year"
la var r_total_der_fica "Total Earnings, FICA jobs, 1st year"
la var r_total_der_nonfica "Total Earnings, non-FICA jobs, 1st year"
la var db_pension "Defined Benefit Pension Plan"
la var dc_pension "Defined Contrib. Pension Plan"
la var pension_in_scope_empl "In Scope for Pension"
la var life_ins_1 "Has Life Insurance"
la var r_fsamt "Food Stamps Amt"
la var hicov "Has Health Insurance"
la var hiemp "Has HI from employer"
la var r_afdcamt "AFDC Amt"

*add afdc indicator
gen afdc_ind = (afdcamt>0) // =1 if positive, 0 o/w
replace afdc_ind = . if mi(afdcamt) // if missing (never happens in SSB implicate), then change to missing here
la var afdc_ind "=1 if receives AFDC"

local dems "sipp_beg_age age_elig male white black otherrace hisp nohs hsdeg somecoll colldeg graddeg married nkids fs_ind no_earn afdc_ind"

******* Section 4: Generate Summary Stats Table*******

svyset person_id [pweight=initwgt]

foreach var in `dems' `lhs' `real_lhs' {

foreach yr in 1996 2001 2004 2008 {

*I'm counting how many households are at the lower, upper and family max bend points and how many total DI households there are in each year.
	*weighted US pop
	*lower bend point
	qui su initwgt if l_`yr' == 1 & year == `yr'
	local lN_`yr'_1=`r(sum)'
	
	*upper bend point
	qui su initwgt if u_`yr' == 1 & year == `yr'
	local uN_`yr'_1=`r(sum)'
	
	*family bend point
	qui su initwgt if f_`yr' == 1 & year == `yr'
	local fN_`yr'_1=`r(sum)'
	
	*DI main sample
	qui su initwgt if year == `yr' & ssdi_no_ssi==1
	local nSDIN_`yr'_1=`r(sum)'
	
	*all non DI
	qui su initwgt if year == `yr' & all_no_ssdi==1
	local nDIN_`yr'_1=`r(sum)'
	
	*full sample
	qui su initwgt if year == `yr'
	local N_`yr'_1=`r(sum)'
	
	*all DI recipients
	qui su initwgt if year == `yr' & all_ssdi==1
	local DIN_`yr'_1=`r(sum)'
	
	*non DI, removing SSI recipients
	qui su initwgt if year == `yr' & no_ssi_no_ssdi==1
	local nSnDIN_`yr'_1=`r(sum)'
	
	*non DI, SSI recipients only
	qui su initwgt if year == `yr' & ssi_no_ssdi==1
	local SnDIN_`yr'_1=`r(sum)'
	
	*serial beneficiaries (SSI and SSDI)
	qui su initwgt if year == `yr' & ssi_and_ssdi==1
	local SDIN_`yr'_1=`r(sum)'
	
	*all SSI beneficiaries
	qui su initwgt if year == `yr' & all_ssi==1
	local SN_`yr'_1=`r(sum)'
	
	*weighted SSB
	
	*first create scale factor to scale down to SSB weighted sum
	*SSB total
	qui count if year == `yr'
	loc ssb_total_`yr' = `r(N)' 
	
	*US pop total
	qui su initwgt if year == `yr'
	loc us_pop_total_`yr'=`r(sum)'
	
	*scale is ratio of the two
	loc scale= `ssb_total_`yr''/`us_pop_total_`yr''
	
	*scale down
	qui su initwgt if l_`yr' == 1 & year == `yr'
	local lN_`yr'_2=`r(sum)'*`scale'
	
	qui su initwgt if u_`yr' == 1 & year == `yr'
	local uN_`yr'_2=`r(sum)'*`scale'
	
	qui su initwgt if f_`yr' == 1 & year == `yr'
	local fN_`yr'_2=`r(sum)'*`scale'
	
	qui su initwgt if year == `yr' & ssdi_no_ssi==1
	local nSDIN_`yr'_2=`r(sum)'*`scale'

	qui su initwgt if year == `yr' & all_no_ssdi==1
	local nDIN_`yr'_2=`r(sum)'*`scale'

	qui su initwgt if year == `yr'
	local N_`yr'_2=`r(sum)'*`scale'
	
	qui su initwgt if year == `yr' & all_ssdi==1
	local DIN_`yr'_2=`r(sum)'*`scale'
	
	qui su initwgt if year == `yr' & no_ssi_no_ssdi==1
	local nSnDIN_`yr'_2=`r(sum)'*`scale'
	
	qui su initwgt if year == `yr' & ssi_no_ssdi==1
	local SnDIN_`yr'_2=`r(sum)'*`scale'
	
	qui su initwgt if year == `yr' & ssi_and_ssdi==1
	local SDIN_`yr'_2=`r(sum)'*`scale'
	
	qui su initwgt if year == `yr' & all_ssi==1
	local SN_`yr'_2=`r(sum)'*`scale'
	
	*unweighted SSB		
	qui count if l_`yr' == 1 & year == `yr'
	local lN_`yr'_3=`r(N)'
	
	qui count if u_`yr' == 1 & year == `yr'
	local uN_`yr'_3=`r(N)'
	
	qui count if f_`yr' == 1 & year == `yr'
	local fN_`yr'_3=`r(N)'
	
	qui count if year == `yr' & ssdi_no_ssi==1
	local nSDIN_`yr'_3=`r(N)'

	qui count if year == `yr' & all_no_ssdi==1
	local nDIN_`yr'_3=`r(N)'

	qui count if year == `yr'
	local N_`yr'_3=`r(N)'
	
	qui count if year == `yr' & all_ssdi==1
	local DIN_`yr'_3=`r(N)'
	
	qui count if year == `yr' & no_ssi_no_ssdi==1
	local nSnDIN_`yr'_3=`r(N)'
	
	qui count if year == `yr' & ssi_no_ssdi==1
	local SnDIN_`yr'_3=`r(N)'
	
	qui count if year == `yr' & ssi_and_ssdi==1
	local SDIN_`yr'_3=`r(N)'
	
	qui count if year == `yr' & all_ssi==1
	local SN_`yr'_3=`r(N)'
}	
* compute means at group level

	*bp1
	qui svy: mean `var', over(group_bp1)
	qui estat sd
	qui mat mean = r(mean)
	qui mat sd = r(sd)
	local lm_`var'=(mean[1,1])
	local lsd_`var'=(sd[1,1])
	quietly _pctile `var' if group_bp1 == 1 [pw=initwgt]
	local lmd_`var'=(r(r1))
	
	*bp2
	quietly svy: mean `var', over(group_bp2)
	quietly estat sd
	qui mat mean = r(mean)
	qui mat sd = r(sd)
	local um_`var'=(mean[1,1])
	local usd_`var'=(sd[1,1])
	quietly _pctile `var' if group_bp2==1 [pw=initwgt]
	local umd_`var'=(r(r1))
	
	*bpf
	quietly svy: mean `var', over(group_bpf)
	quietly estat sd
	qui mat mean = r(mean)
	qui mat sd = r(sd)
	local fm_`var'=(mean[1,1])
	local fsd_`var'=(sd[1,1])
	quietly _pctile `var' if group_bpf==1 [pw=initwgt]
	local fmd_`var'=(r(r1))
	
	*DI main sample
	quietly svy: mean `var', over(ssdi_no_ssi)
	quietly estat sd
	qui mat mean = r(mean)
	qui mat sd = r(sd)
	local nSDIm_`var'=(mean[1,1])
	local nSDIsd_`var'=(sd[1,1])
	quietly _pctile `var' [pw=initwgt]
	local nSDImd_`var'=(r(r1))
	
	*Non DI
	quietly svy: mean `var', over(all_no_ssdi)
	quietly estat sd
	qui mat mean = r(mean)
	qui mat sd = r(sd)
	local nDIm_`var'=(mean[1,1])
	local nDIsd_`var'=(sd[1,1])
	quietly _pctile `var' [pw=initwgt]
	local nDImd_`var'=(r(r1))

	*Full Sample
	quietly svy: mean `var'
	quietly estat sd
	qui mat mean = r(mean)
	qui mat sd = r(sd)
	local m_`var'=(mean[1,1])
	local sd_`var'=(sd[1,1])
	quietly _pctile `var' [pw=initwgt]
	local md_`var'=(r(r1))
	
	*all DI
	quietly svy: mean `var', over(all_ssdi)
	quietly estat sd
	qui mat mean = r(mean)
	qui mat sd = r(sd)
	local DIm_`var'=(mean[1,1])
	local DIsd_`var'=(sd[1,1])
	quietly _pctile `var' [pw=initwgt]
	local DImd_`var'=(r(r1))
	
	*Non DI, removing SSI recipients
	quietly svy: mean `var', over(no_ssi_no_ssdi)
	quietly estat sd
	qui mat mean = r(mean)
	qui mat sd = r(sd)
	local nSnDIm_`var'=(mean[1,1])
	local nSnDIsd_`var'=(sd[1,1])
	quietly _pctile `var' [pw=initwgt]
	local nSnDImd_`var'=(r(r1))
	
	*Non DI, SSI only
	quietly svy: mean `var', over(ssi_no_ssdi)
	quietly estat sd
	qui mat mean = r(mean)
	qui mat sd = r(sd)
	local SnDIm_`var'=(mean[1,1])
	local SnDIsd_`var'=(sd[1,1])
	quietly _pctile `var' [pw=initwgt]
	local SnDImd_`var'=(r(r1))
	
	*SSI and SSDI
	quietly svy: mean `var', over(ssi_and_ssdi)
	quietly estat sd
	qui mat mean = r(mean)
	qui mat sd = r(sd)
	local SDIm_`var'=(mean[1,1])
	local SDIsd_`var'=(sd[1,1])
	quietly _pctile `var' [pw=initwgt]
	local SDImd_`var'=(r(r1))
	
	*all SSI
	quietly svy: mean `var', over(all_ssi)
	quietly estat sd
	qui mat mean = r(mean)
	qui mat sd = r(sd)
	local Sm_`var'=(mean[1,1])
	local Ssd_`var'=(sd[1,1])
	quietly _pctile `var' [pw=initwgt]
	local Smd_`var'=(r(r1))
	
}

/* And finally, I'm just going to cycle through all of my variables and for each one, average the summary statistics across the 4 surveys. 
And I'll place the average in the cell corresponding to that variable and that group. */
qui putexcel set sumstats_SSDI_SSB_${gsnum}${impflag}, mod

	qui putexcel A1 = ("Variable"), bold underline fpat(solid, white)
	qui putexcel B1 = ("Variable Label"), bold underline fpat(solid, white)
	qui putexcel D1 = ("Lower Bend Point"), hcenter bold underline fpat(solid, white)
	qui putexcel H1 = ("Upper Bend Point"), hcenter bold underline fpat(solid, white)
	qui putexcel L1 = ("Family Max Bend Point"), hcenter bold underline fpat(solid, white)
	qui putexcel P1 = ("DI Sample"), hcenter bold underline fpat(solid, white)
	qui putexcel T1 = ("All Non-DI Recipients"), hcenter bold underline fpat(solid, white)
	qui putexcel X1 = ("Full Sample"), hcenter bold underline fpat(solid, white)
	qui putexcel AB1 = ("All DI"), hcenter bold underline fpat(solid, white)
	qui putexcel AF1 = ("Non-DI, No SSI"), hcenter bold underline fpat(solid, white)
	qui putexcel AJ1 = ("Non-DI, SSI Only"), hcenter bold underline fpat(solid, white)
	qui putexcel AN1 = ("SSI-SSDI Dual Beneficiaries"), hcenter bold underline fpat(solid, white)
	qui putexcel AR1 = ("All SSI"), hcenter bold underline fpat(solid, white)

	qui putexcel D3 = ("Median"), italic right fpat(solid, white)
	qui putexcel E3 = ("Mean"), italic right fpat(solid, white)
	qui putexcel F3 = ("St Dev"), italic right fpat(solid, white)

	qui putexcel H3 = ("Median"), italic right fpat(solid, white)
	qui putexcel I3 = ("Mean"), italic right fpat(solid, white)
	qui putexcel J3 = ("St Dev"), italic right fpat(solid, white)
	
	qui putexcel L3 = ("Median"), italic right fpat(solid, white)
	qui putexcel M3 = ("Mean"), italic right fpat(solid, white)
	qui putexcel N3 = ("St Dev"), italic right fpat(solid, white)

	qui putexcel P3 = ("Median"), italic right fpat(solid, white)
	qui putexcel Q3 = ("Mean"), italic right fpat(solid, white)
	qui putexcel R3 = ("St Dev"), italic right fpat(solid, white)
	
	qui putexcel T3 = ("Median"), italic right fpat(solid, white)
	qui putexcel U3 = ("Mean"), italic right fpat(solid, white)
	qui putexcel V3 = ("St Dev"), italic right fpat(solid, white)

	qui putexcel X3 = ("Median"), italic right fpat(solid, white)
	qui putexcel Y3 = ("Mean"), italic right fpat(solid, white)
	qui putexcel Z3 = ("St Dev"), italic right fpat(solid, white)
	
	qui putexcel AB3 = ("Median"), italic right fpat(solid, white)
	qui putexcel AC3 = ("Mean"), italic right fpat(solid, white)
	qui putexcel AD3 = ("St Dev"), italic right fpat(solid, white)
	
	qui putexcel AF3 = ("Median"), italic right fpat(solid, white)
	qui putexcel AG3 = ("Mean"), italic right fpat(solid, white)
	qui putexcel AH3 = ("St Dev"), italic right fpat(solid, white)
	
	qui putexcel AJ3 = ("Median"), italic right fpat(solid, white)
	qui putexcel AK3 = ("Mean"), italic right fpat(solid, white)
	qui putexcel AL3 = ("St Dev"), italic right fpat(solid, white)
	
	qui putexcel AN3 = ("Median"), italic right fpat(solid, white)
	qui putexcel AO3 = ("Mean"), italic right fpat(solid, white)
	qui putexcel AP3 = ("St Dev"), italic right fpat(solid, white)
	
	qui putexcel AR3 = ("Median"), italic right fpat(solid, white)
	qui putexcel AS3 = ("Mean"), italic right fpat(solid, white)
	qui putexcel AT3 = ("St Dev"), italic right fpat(solid, white)
	
	qui putexcel D2 = ("N="), right fpat(solid, white)
	qui putexcel H2 = ("N="), right fpat(solid, white)
	qui putexcel L2 = ("N="), right fpat(solid, white)
	qui putexcel P2 = ("N="), right fpat(solid, white)
	qui putexcel T2 = ("N="), right fpat(solid, white)
	qui putexcel X2 = ("N="), right fpat(solid, white)
	qui putexcel AB2 = ("N="), right fpat(solid, white)
	qui putexcel AF2 = ("N="), right fpat(solid, white)
	qui putexcel AJ2 = ("N="), right fpat(solid, white)
	qui putexcel AN2 = ("N="), right fpat(solid, white)
	qui putexcel AR2 = ("N="), right fpat(solid, white)

	*US pop level
	qui putexcel E2 = (((`lN_1996_1'+`lN_2001_1'+`lN_2004_1'+`lN_2008_1'))), left fpat(solid, white) nformat("#.#")
	qui putexcel I2 = (((`uN_1996_1'+`uN_2001_1'+`uN_2004_1'+`uN_2008_1'))), left fpat(solid, white) nformat("#.#")
	qui putexcel M2 = (((`fN_1996_1'+`fN_2001_1'+`fN_2004_1'+`fN_2008_1'))), left fpat(solid, white) nformat("#.#")
	qui putexcel Q2 = (((`nSDIN_1996_1'+`nSDIN_2001_1'+`nSDIN_2004_1'+`nSDIN_2008_1'))), left fpat(solid, white) nformat("#.#")
	qui putexcel U2 = (((`nDIN_1996_1'+`nDIN_2001_1'+`nDIN_2004_1'+`nDIN_2008_1'))), left fpat(solid, white) nformat("#.#")
	qui putexcel Y2 = (((`N_1996_1'+`N_2001_1'+`N_2004_1'+`N_2008_1'))), left fpat(solid, white) nformat("#.#")
	qui putexcel AC2 = (((`DIN_1996_1'+`DIN_2001_1'+`DIN_2004_1'+`DIN_2008_1'))), left fpat(solid, white) nformat("#.#")
	qui putexcel AG2 = (((`nSnDIN_1996_1'+`nSnDIN_2001_1'+`nSnDIN_2004_1'+`nSnDIN_2008_1'))), left fpat(solid, white) nformat("#.#")
	qui putexcel AK2 = (((`SnDIN_1996_1'+`SnDIN_2001_1'+`SnDIN_2004_1'+`SnDIN_2008_1'))), left fpat(solid, white) nformat("#.#")
	qui putexcel AO2 = (((`SDIN_1996_1'+`SDIN_2001_1'+`SDIN_2004_1'+`SDIN_2008_1'))), left fpat(solid, white) nformat("#.#")
	qui putexcel AS2 = (((`SN_1996_1'+`SN_2001_1'+`SN_2004_1'+`SN_2008_1'))), left fpat(solid, white) nformat("#.#")
	
	*SSB weighted
	qui putexcel F2 = (((`lN_1996_2'+`lN_2001_2'+`lN_2004_2'+`lN_2008_2'))), left fpat(solid, white) nformat("#.#")
	qui putexcel J2 = (((`uN_1996_2'+`uN_2001_2'+`uN_2004_2'+`uN_2008_2'))), left fpat(solid, white) nformat("#.#")
	qui putexcel N2 = (((`fN_1996_2'+`fN_2001_2'+`fN_2004_2'+`fN_2008_2'))), left fpat(solid, white) nformat("#.#")
	qui putexcel R2 = (((`nSDIN_1996_2'+`nSDIN_2001_2'+`nSDIN_2004_2'+`nSDIN_2008_2'))), left fpat(solid, white) nformat("#.#")
	qui putexcel V2 = (((`nDIN_1996_2'+`nDIN_2001_2'+`nDIN_2004_2'+`nDIN_2008_2'))), left fpat(solid, white) nformat("#.#")
	qui putexcel Z2 = (((`N_1996_2'+`N_2001_2'+`N_2004_2'+`N_2008_2'))), left fpat(solid, white) nformat("#.#")
	qui putexcel AD2 = (((`DIN_1996_2'+`DIN_2001_2'+`DIN_2004_2'+`DIN_2008_2'))), left fpat(solid, white) nformat("#.#")
	qui putexcel AH2 = (((`nSnDIN_1996_2'+`nSnDIN_2001_2'+`nSnDIN_2004_2'+`nSnDIN_2008_2'))), left fpat(solid, white) nformat("#.#")
	qui putexcel AL2 = (((`SnDIN_1996_2'+`SnDIN_2001_2'+`SnDIN_2004_2'+`SnDIN_2008_2'))), left fpat(solid, white) nformat("#.#")
	qui putexcel AP2 = (((`SDIN_1996_2'+`SDIN_2001_2'+`SDIN_2004_2'+`SDIN_2008_2'))), left fpat(solid, white) nformat("#.#")
	qui putexcel AT2 = (((`SN_1996_2'+`SN_2001_2'+`SN_2004_2'+`SN_2008_2'))), left fpat(solid, white) nformat("#.#")
		
	*SSB unweighted
	qui putexcel G2 = (((`lN_1996_3'+`lN_2001_3'+`lN_2004_3'+`lN_2008_3'))), left fpat(solid, white) nformat("#.#")
	qui putexcel K2 = (((`uN_1996_3'+`uN_2001_3'+`uN_2004_3'+`uN_2008_3'))), left fpat(solid, white) nformat("#.#")
	qui putexcel O2 = (((`fN_1996_3'+`fN_2001_3'+`fN_2004_3'+`fN_2008_3'))), left fpat(solid, white) nformat("#.#")
	qui putexcel S2 = (((`nSDIN_1996_3'+`nSDIN_2001_3'+`nSDIN_2004_3'+`nSDIN_2008_3'))), left fpat(solid, white) nformat("#.#")
	qui putexcel W2 = (((`nDIN_1996_3'+`nDIN_2001_3'+`nDIN_2004_3'+`nDIN_2008_3'))), left fpat(solid, white) nformat("#.#")
	qui putexcel AA2 = (((`N_1996_3'+`N_2001_3'+`N_2004_3'+`N_2008_3'))), left fpat(solid, white) nformat("#.#")
	qui putexcel AE2 = (((`DIN_1996_3'+`DIN_2001_3'+`DIN_2004_3'+`DIN_2008_3'))), left fpat(solid, white) nformat("#.#")
	qui putexcel AI2 = (((`nSnDIN_1996_3'+`nSnDIN_2001_3'+`nSnDIN_2004_3'+`nSnDIN_2008_3'))), left fpat(solid, white) nformat("#.#")
	qui putexcel AM2 = (((`SnDIN_1996_3'+`SnDIN_2001_3'+`SnDIN_2004_3'+`SnDIN_2008_3'))), left fpat(solid, white) nformat("#.#")
	qui putexcel AQ2 = (((`SDIN_1996_3'+`SDIN_2001_3'+`SDIN_2004_3'+`SDIN_2008_3'))), left fpat(solid, white) nformat("#.#")
	qui putexcel AU2 = (((`SN_1996_3'+`SN_2001_3'+`SN_2004_3'+`SN_2008_3'))), left fpat(solid, white) nformat("#.#")

local row = 4

foreach var in `dems' `lhs' `real_lhs' {

	qui putexcel A`row' = ("`var'")
	local varlabel : var label `var'
	qui putexcel B`row' = ("`varlabel'")
	
	qui putexcel D`row' = (`lmd_`var''), fpat(solid, white) nformat(number_d2)
	qui putexcel E`row' = (`lm_`var''), fpat(solid, white) nformat(number_d2)
	qui putexcel F`row' = (`lsd_`var''), fpat(solid, white) nformat(number_d2)
	
	qui putexcel H`row' = (`umd_`var''), fpat(solid, white) nformat(number_d2)
	qui putexcel I`row' = (`um_`var''), fpat(solid, white) nformat(number_d2)
	qui putexcel J`row' = (`usd_`var''), fpat(solid, white) nformat(number_d2)

	qui putexcel L`row' = (`fmd_`var''), fpat(solid, white) nformat(number_d2)
	qui putexcel M`row' = (`fm_`var''), fpat(solid, white) nformat(number_d2)
	qui putexcel N`row' = (`fsd_`var''), fpat(solid, white) nformat(number_d2)
	
	qui putexcel P`row' = (`nSDImd_`var''), fpat(solid, white) nformat(number_d2)
	qui putexcel Q`row' = (`nSDIm_`var''), fpat(solid, white) nformat(number_d2)
	qui putexcel R`row' = (`nSDIsd_`var''), fpat(solid, white) nformat(number_d2)

	qui putexcel T`row' = (`nDImd_`var''), fpat(solid, white) nformat(number_d2)
	qui putexcel U`row' = (`nDIm_`var''), fpat(solid, white) nformat(number_d2)
	qui putexcel V`row' = (`nDIsd_`var''), fpat(solid, white) nformat(number_d2)

	qui putexcel X`row' = (`md_`var''), fpat(solid, white) nformat(number_d2)
	qui putexcel Y`row' = (`m_`var''), fpat(solid, white) nformat(number_d2)
	qui putexcel Z`row' = (`sd_`var''), fpat(solid, white) nformat(number_d2)

	qui putexcel AB`row' = (`DImd_`var''), fpat(solid, white) nformat(number_d2)
	qui putexcel AC`row' = (`DIm_`var''), fpat(solid, white) nformat(number_d2)
	qui putexcel AD`row' = (`DIsd_`var''), fpat(solid, white) nformat(number_d2)

	qui putexcel AF`row' = (`nSnDImd_`var''), fpat(solid, white) nformat(number_d2)
	qui putexcel AG`row' = (`nSnDIm_`var''), fpat(solid, white) nformat(number_d2)
	qui putexcel AH`row' = (`nSnDIsd_`var''), fpat(solid, white) nformat(number_d2)

	qui putexcel AJ`row' = (`SnDImd_`var''), fpat(solid, white) nformat(number_d2)
	qui putexcel AK`row' = (`SnDIm_`var''), fpat(solid, white) nformat(number_d2)
	qui putexcel AL`row' = (`SnDIsd_`var''), fpat(solid, white) nformat(number_d2)

	qui putexcel AN`row' = (`SDImd_`var''), fpat(solid, white) nformat(number_d2)
	qui putexcel AO`row' = (`SDIm_`var''), fpat(solid, white) nformat(number_d2)
	qui putexcel AP`row' = (`SDIsd_`var''), fpat(solid, white) nformat(number_d2)

	qui putexcel AR`row' = (`Smd_`var''), fpat(solid, white) nformat(number_d2)
	qui putexcel AS`row' = (`Sm_`var''), fpat(solid, white) nformat(number_d2)
	qui putexcel AT`row' = (`Ssd_`var''), fpat(solid, white) nformat(number_d2)

	local ++ row
}
}

*SSDI only variables

loc ssionly reason_mental reason_musculo reason_other inc_allDI

foreach var in `ssionly' {

* compute means at group level

	*bp1
	qui svy: mean `var', over(group_bp1)
	qui estat sd
	qui mat mean = r(mean)
	qui mat sd = r(sd)
	local lm_`var'=(mean[1,1])
	local lsd_`var'=(sd[1,1])
	quietly _pctile `var' if group_bp1 == 1 [pw=initwgt]
	local lmd_`var'=(r(r1))
	
	*bp2
	quietly svy: mean `var', over(group_bp2)
	quietly estat sd
	qui mat mean = r(mean)
	qui mat sd = r(sd)
	local um_`var'=(mean[1,1])
	local usd_`var'=(sd[1,1])
	quietly _pctile `var' if group_bp2==1 [pw=initwgt]
	local umd_`var'=(r(r1))
	
	*bpf
	quietly svy: mean `var', over(group_bpf)
	quietly estat sd
	qui mat mean = r(mean)
	qui mat sd = r(sd)
	local fm_`var'=(mean[1,1])
	local fsd_`var'=(sd[1,1])
	quietly _pctile `var' if group_bpf==1 [pw=initwgt]
	local fmd_`var'=(r(r1))
	
	*DI main sample
	quietly svy: mean `var', over(ssdi_no_ssi)
	quietly estat sd
	qui mat mean = r(mean)
	qui mat sd = r(sd)
	local nSDIm_`var'=(mean[1,1])
	local nSDIsd_`var'=(sd[1,1])
	quietly _pctile `var' [pw=initwgt]
	local nSDImd_`var'=(r(r1))
	
	*Full Sample
	quietly svy: mean `var'
	quietly estat sd
	qui mat mean = r(mean)
	qui mat sd = r(sd)
	local m_`var'=(mean[1,1])
	local sd_`var'=(sd[1,1])
	quietly _pctile `var' [pw=initwgt]
	local md_`var'=(r(r1))
	
	*all DI
	quietly svy: mean `var', over(all_ssdi)
	quietly estat sd
	qui mat mean = r(mean)
	qui mat sd = r(sd)
	local DIm_`var'=(mean[1,1])
	local DIsd_`var'=(sd[1,1])
	quietly _pctile `var' [pw=initwgt]
	local DImd_`var'=(r(r1))
		
	*SSI and SSDI
	quietly svy: mean `var', over(ssi_and_ssdi)
	quietly estat sd
	qui mat mean = r(mean)
	qui mat sd = r(sd)
	local SDIm_`var'=(mean[1,1])
	local SDIsd_`var'=(sd[1,1])
	quietly _pctile `var' [pw=initwgt]
	local SDImd_`var'=(r(r1))
	
	*all SSI
	quietly svy: mean `var', over(all_ssi)
	quietly estat sd
	qui mat mean = r(mean)
	qui mat sd = r(sd)
	local Sm_`var'=(mean[1,1])
	local Ssd_`var'=(sd[1,1])
	quietly _pctile `var' [pw=initwgt]
	local Smd_`var'=(r(r1))
	
}

di `row'

foreach var in `ssionly' {

	qui putexcel A`row' = ("`var'")
	local varlabel : var label `var'
	qui putexcel B`row' = ("`varlabel'")
	
	qui putexcel D`row' = (`lmd_`var''), fpat(solid, white) nformat(number_d2)
	qui putexcel E`row' = (`lm_`var''), fpat(solid, white) nformat(number_d2)
	qui putexcel F`row' = (`lsd_`var''), fpat(solid, white) nformat(number_d2)
	
	qui putexcel H`row' = (`umd_`var''), fpat(solid, white) nformat(number_d2)
	qui putexcel I`row' = (`um_`var''), fpat(solid, white) nformat(number_d2)
	qui putexcel J`row' = (`usd_`var''), fpat(solid, white) nformat(number_d2)

	qui putexcel L`row' = (`fmd_`var''), fpat(solid, white) nformat(number_d2)
	qui putexcel M`row' = (`fm_`var''), fpat(solid, white) nformat(number_d2)
	qui putexcel N`row' = (`fsd_`var''), fpat(solid, white) nformat(number_d2)
	
	qui putexcel P`row' = (`nSDImd_`var''), fpat(solid, white) nformat(number_d2)
	qui putexcel Q`row' = (`nSDIm_`var''), fpat(solid, white) nformat(number_d2)
	qui putexcel R`row' = (`nSDIsd_`var''), fpat(solid, white) nformat(number_d2)

	qui putexcel T`row' = (.), fpat(solid, white) nformat(number_d2)
	qui putexcel U`row' = (.), fpat(solid, white) nformat(number_d2)
	qui putexcel V`row' = (.), fpat(solid, white) nformat(number_d2)

	qui putexcel X`row' = (`md_`var''), fpat(solid, white) nformat(number_d2)
	qui putexcel Y`row' = (`m_`var''), fpat(solid, white) nformat(number_d2)
	qui putexcel Z`row' = (`sd_`var''), fpat(solid, white) nformat(number_d2)

	qui putexcel AB`row' = (`DImd_`var''), fpat(solid, white) nformat(number_d2)
	qui putexcel AC`row' = (`DIm_`var''), fpat(solid, white) nformat(number_d2)
	qui putexcel AD`row' = (`DIsd_`var''), fpat(solid, white) nformat(number_d2)

	qui putexcel AF`row' = (.), fpat(solid, white) nformat(number_d2)
	qui putexcel AG`row' = (.), fpat(solid, white) nformat(number_d2)
	qui putexcel AH`row' = (.), fpat(solid, white) nformat(number_d2)

	qui putexcel AJ`row' = (.), fpat(solid, white) nformat(number_d2)
	qui putexcel AK`row' = (.), fpat(solid, white) nformat(number_d2)
	qui putexcel AL`row' = (.), fpat(solid, white) nformat(number_d2)

	qui putexcel AN`row' = (`SDImd_`var''), fpat(solid, white) nformat(number_d2)
	qui putexcel AO`row' = (`SDIm_`var''), fpat(solid, white) nformat(number_d2)
	qui putexcel AP`row' = (`SDIsd_`var''), fpat(solid, white) nformat(number_d2)

	qui putexcel AR`row' = (`Smd_`var''), fpat(solid, white) nformat(number_d2)
	qui putexcel AS`row' = (`Sm_`var''), fpat(solid, white) nformat(number_d2)
	qui putexcel AT`row' = (`Ssd_`var''), fpat(solid, white) nformat(number_d2)

	local ++ row
}

erase famids.dta

*now we show which cells use less than 10 observations by cross-tabulating the dummy variables across each sample:
loc indicators married reason_mental reason_musculo reason_other fs_ind no_earn inc_allDI male white black otherrace hisp nohs hsdeg somecoll ///
	colldeg graddeg own_home db_pension dc_pension pension_in_scope_empl life_ins_1 hicov hiemp afdc_ind
	
foreach var in `indicators' {
di "Variable: `var'"
di "Full Sample"
tab `var'
di "Lower Bend Pt"
tab `var' if group_bp1==1
di "Upper Bend Pt"
tab `var' if group_bp2==1
di "Family Max Bend Pt"
tab `var' if group_bpf==1
di "All SSDI"
tab `var' if all_ssdi == 1
di "All No SSDI"
tab `var' if all_no_ssdi == 1
di "All SSDI without SSI"
tab `var' if ssdi_no_ssi == 1
di "All SSI without SSDI"
tab `var' if ssi_no_ssdi == 1
di "All SSI with SSDI"
tab `var' if ssi_and_ssdi == 1
di "All Non-SSI without SSDI"
tab `var' if no_ssi_no_ssdi == 1
di "All SSI"
tab `var' if all_ssi == 1
di "							"
}
